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Recently the general form of a translation-covariant quantum Boltzmann equation has been de- 
rived which describes the dynamics of a tracer particle in a quantum gas. We develop a stochastic 
wave function algorithm that enables full three-dimensional Monte Carlo simulations of this equa- 
tion. The simulation method is used to study the approach to equilibrium for various scattering 
cross sections and to determine dynamical deviations from Gaussian statistics through an investi- 
gation of higher-order cumulants. Moreover, we examine the loss of coherence of superpositions of 
momentum eigenstates and determine the corresponding decoherence time scales to quantify the 
transition from quantum to classical behavior of the state of the test particle. 

PACS numbers: 03.65.Yz, 02.70.Ss, 05.20.Dd, 47.45.Ab 



^ ■ 

-f— > 



> 

CO 

in 

(N 

o 
o 



X 



I. INTRODUCTION 

In recent times major efforts have been devoted to the 
study and understanding of the dynamics of open systems 
[l[, both in order to give a realistic, quantitative descrip- 
tion of the time evolution of a quantum system coupled 
to a generally larger system considered as environment, 
as well as with the aim to engineer suitable environments 
driving the dynamics of the quantum system according 
to the will of the experimenter. When considering an 
open system one can either use a phenomenological de- 
scription for the system-environment interaction, as well 
as for the environment itself, or rely on a strictly mi- 
croscopic description of the physical system considered. 
While the first approach can be more viable and flexible, 
the second is clearly of more fundamental nature. In the 
present paper we will focus on the second approach, per- 
forming a numerical study of a recently observed @, HJ 
quantum master equation for the description of the mo- 
tion of a quantum test particle in a gas. Such a master 
equation is the quantum version of the classical linear 
Boltzmann equation and gives a microscopic description 
of the dynamics of the test particle, only relying on the 
gas properties and on the detailed expression of the in- 
teraction between test particle and gas particles. The 
quantum linear Boltzmann equation is characterized by 
its covariance under translations and has been used 
in a simplified form for the quantitative description of 
experiments on collisional decoherence [1, H, 0] ■ Quanti- 
tative experiments, testing the transition from the quan- 
tum to the classical world, are in fact typical situations 
in which a truly microscopic description of the physical 
system of interest is mandatory. 

The translation-covariant quantum linear Boltzmann 



equation that will be studied here may be written in 
the form of a quantum master equation for the time- 
dependent density matrix p(t) of the test particle, 



dt 



p(t) = Cp{t). 



(1) 
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The infinitesimal generator C represents a superopera- 
tor in Lindblad form [1, Q . It is well known that such 
a quantum master equation with Lindblad structure al- 
lows an unravelling through a stochastic process for the 
state vector in the particle's Hilbert space. Here, we 
will concentrate on the so-called quantum jump method 
in which the state vector follows a piecewise determin- 
istic process, consisting of smooth, deterministic evolu- 
tion periods and discontinuous, random quantum jumps 
[Tol [ill H3 . [l3| . It will be demonstrated that the applica- 
tion of this method to the quantum Boltzmann equation 
(fTj) is indeed feasible and leads to a simple and numeri- 
cally efficient three-dimensional Monte Carlo simulation 
technique for the dynamical behavior of the particle. 

By means of this technique we will study in particu- 
lar relaxation properties of this master equation for var- 
ious microscopic scattering cross sections, as well as de- 
viations from Gaussian statistics. Given the quantum 
nature of the equation we will also investigate the time 
evolution of quantum superposition states. 

The paper is organized as follows. In Sec. [Til we briefly 
introduce the quantum linear Boltzmann equation to- 
gether with its basic features, while in Sec. Mil we show 
how to apply the Monte Carlo wave function method to 
this master equation, which has a very intuitive physi- 
cal meaning when working in the momentum represen- 
tation. In Sec. IIVI we describe the numerical algorithms 
used and present the simulation results. Besides studying 
relaxation dynamics and deviation from Gaussian statis- 
tics, we will focus on decoherence effects in momentum 
space, comparing relaxation and decoherence rates and 
providing analytical estimates for the latter. Finally, in 
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Sec.|V]we comment on our results and point out possible 
extensions and generalizations of the work. 



II. QUANTUM MASTER EQUATION 



Q and energy transfer J5(Q, P) of the time dependent 
density-density autocorrelation function of the medium, 

S(Q,E) = ^-JdtJ dXei^-Q^) G{X:t)i (4) 



We briefly introduce the quantum master equation 
that we want to study numerically, together with its most 
relevant features, referring the reader to 0, [H, 0, [HL [H| 
for a more detailed presentation. The explicit form of the 
Lindblad generator of Eq. fl} is given by 



Cp(t) = --[H ,p(t)} 



dQcr(Q) 



-{5(Q,P),p(t)} 



(2) 



where Hq = P 2 /2M is the kinetic energy of the test 
particle with mass M, X and P are its position and mo- 
mentum operators respectively, n gas is the number den- 
sity of the gas, m is the mass of the gas particles, and 
m* = mM/(m + M) denotes the reduced mass. 

The master equation fl]) with the generator ([2]) is the 
quantum version of the classical linear Boltzmann equa- 
tion, provided collisions can be described in Born approx- 
imation, or the scattering cross section only depends on 
the momentum transfer experienced by the test particle 
in the collision. We denote such a scattering cross sec- 
tion by c(Q). A more general expression of the quantum 
linear Boltzmann equation, including the full scattering 
cross section on general grounds, has been obtained in 
3, and relies on the appearance of the operator- valued 
scattering amplitude. It coincides with Eq. ([T]) in Born 
approximation or when the full scattering cross section 
only depends on the momentum transfer. Such a master 
equation describes the motion of a quantum test par- 
ticle interacting through collisions with a dilute gas of 
environmental particles. The positive quantity i5(Q,P), 
here appearing operator-valued, is a two-point correla- 
tion function of the gas related to the spectrum of its 
density fluctuations. It is most often expressed in its de- 
pendence on the momentum transfer Q and on the energy 
transfer E(Q, P) = Q 2 /2M + P • Q/M characterizing a 
collision in which the test particle gains a momentum Q, 
changing its momentum from P to P + Q. For the case of 
a free gas of particles described by a Maxwell-Boltzmann 
distribution it is explicitly given by 



S(Q,P) = S(Q,£(Q,P)) 



' f3m 1 



(3) 



(3 (2m£:(Q,P) + Q 2 ) 2 
8m Q 2 



where f3 = is the inverse temperature of the gas. 

In full generality this two-point correlation function 
known as dynamic structure factor [TrI [l7j is given by 
the Fourier transform with respect to momentum transfer 



where 

G(X, t) = 1 J dY (N(Y)N(X + Y, t)) (5) 

describes density correlations in the gas. While for the 
general case of an interacting system of particles an exact 
evaluation is obviously not feasible, the dynamic struc- 
ture factor has several important properties that are help- 
ful in the construction of a phenomenological ansatz. 
An important property of the dynamic structure factor 
granting the existence of the expected canonical station- 
ary solution of ([!]) is the so-called detailed balance con- 
dition according to which 



S(Q,E)=e- 0E S(-Q,-E). 



(6) 



Another crucial feature of Eq. ((2]) is its covariance 
under translations. Considering the unitary represen- 
tation t/(a) = exp(— ia ■ P/h), a € R 3 , of the group of 
three-dimensional space translations in the test particle's 
Hilbert space, one has that 



C[U(R)p(t)W(B)] = U(a)£[p(t)]U<(a). 



(7) 



In fact, Eq. $2$ complies with the general structure 
of translation-covariant master equations obtained by 
Holevo [Tq |. providing a physically relevant example of 
this general mathematical structure. The property of co- 
variance reflects the underlying symmetry under transla- 
tions, arising because we are considering a homogeneous 
gas and the interaction potential between test particle 
and gas particles only depends on the relative distance 
between the two. An important consequence of this prop- 
erty, that we shall exploit in the simulations carried out 
in Sec. IIV1 is the fact that the algebra generated by the 
momentum operators P is invariant under the time evo- 
lution described by C. 

The master equation @ is an operator equation, co- 
inciding with the classical linear Boltzmann equation as 
far as the diagonal matrix elements in the momentum 
representation are concerned, but also describing quan- 
tum coherences corresponding to the off-diagonal matrix 
elements, as well as the time evolution of highly non clas- 
sical motional states, such as superposition states. The 
Lindblad structure of such a quantum linear Boltzmann 
equation apart from preservation of trace and positivity 
is of great importance in that it implies the possibility 
to consider suitable stochastic unravellings leading to ef- 
ficient Monte Carlo simulations of the time evolution of 
different quantities of physical interest, despite the high 
complexity of the problem. 
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III. MONTE CARLO WAVE FUNCTION 
METHOD 

We give a short description of the standard stochastic 
jump unravelling [l(| [U [H, G3 of the master equation 
([1]). For further details on the theory and the numer- 
ical implementation see Ref. [l[ and references therein. 
Introducing the Lindblad operators 



L(Q) = e^ x fe<7(Q)S(Q,P) (8) 
V mi 

we can write the Boltzmann equation as follows, 



- P (t) = --[H , P (t)] 



dQ 



L(Q)p(t)Lt(Q) - - {it(Q)i( Q ) )P ( t )} 



(9) 



This equation leads to a stochastic unravelling through 
a piecewise deterministic process in Hilbert space. This 
means that the realizations \ip(t)) of the process consist 
of deterministic evolution periods (deterministic drift) 
which are interrupted by discontinuous changes of the 
state vector (quantum jumps). 



A. Simulation algorithm 

The realizations of the piecewise deterministic process 
are defined through the following algorithm. In between 
two jumps the state vector follows a deterministic time- 
evolution which is given by the nonlinear Schrodinger 
equation 



-I 



with the non-Hcrmitian Hamiltonian 



H cS = H - y J rfQL t (Q)i(Q). 



(10) 



(11) 



Suppose that at time to a jump into some state \ip(to)} = 
\ip) occurred. The total rate for jumps out of this state 
is given by 



dQ (V|£ t (Q)i(Q)|^). 



(12) 



Hence, the next jump will take place at time to+T, where 
t is a stochastic time step which follows the cumulative 
distribution function 



F(t) 



1 



\exp[-iH eS T/h]\ib}\\ 2 



(13) 



This is the waiting time distribution, i. e., F(t) rep- 
resents the probability that the next jump takes place 



somewhere in the interval (to, to + r). Employing the in- 
version method, for instance, one determines the stochas- 
tic time step r by solving the equation 



\ex V [-iH eS T/ti}\^)\\ 2 = V 



(14) 



for t, where r\ is a random number uniformly distributed 
over the interval (0, 1). 

Once the random time step has been determined one 
carries out a jump of the state vector \ip(t + t)) = \tp) 
at time to + t by the replacement 



L(QM) 

wmmw 



(15) 



The momentum transfer Q is to be drawn from the prob- 
ability density 



R(Q) = 

which is normalized as 



MLt(Q)£(Q)|y,) 



J dQ R(Q) = 1. 



(16) 



(17) 



The process thus defined represents a stochastic unrav- 
elling of the quantum master equation in the sense that 
the expectation value 



P (t) = E[m))(m\] 



(18) 



yields a solution of Eq. ([5]). In the Monte Carlo wave 
function method one numerically generates large samples 
of realizations and estimates all desired quantities with 
the help of appropriate sample averages. 



B. Momentum representation 

As already mentioned in Sec. [Til an important property 
of Eq. @ is its covariance under the action of the trans- 
lation group. This implies that the algebra generated 
by the three commuting momentum operators P, i. e. 
the generators of translations, is left invariant under the 
action of the master equation. A function of the momen- 
tum operators goes over with elapsing time to another 
function of the momentum operators only, when evolv- 
ing according to the master equation ([9]) . The considered 
unravelling of Eq. ^ preserves this property in the sense 
that linear combinations of improper eigenvectors |P) of 
the three commuting momentum operators are preserved 
in form in each single realization, which is of great advan- 
tage in the simulations. In fact H c g given by (|11[) is only a 
function of the momentum operators, thus simply acting 
in a multiplicative way on the momentum eigenvectors, 
while the jumps effected by the Lindblad operators L(Q) 
according to (|15|) simply correspond to shifts 

|P)^|P + Q>. 
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Given a master equation covariant under an Abelian 
symmetry group it is generally true that the algebra gen- 
erated by the commuting self-adjoint operators which act 
as generators of the symmetry is left invariant under time 
evolution. Correspondingly one can consider unravcllings 
leaving invariant in form linear combinations of common 
eigenvectors of the generators of the symmetry, where 
the jumps only lead to a shift between different eigenvec- 
tors. For an initial state given by such an eigenvector the 
stochastic unravelling leads to a pure jump process. Con- 
sider for example the master equation for the damped 
harmonic oscillator with Lindblad o per ators a and , 
covariant under the group U(l) [191 l20l. l2ll . [22l |. where 
the generator of the symmetry is the number operator 
N = a^a. In this case one can consider a stochastic un- 
ravelling given by a suitable piecewise deterministic pro- 
cess, where the jumps effected by the Lindblad operators 
a and aJ in the single realizations are simply given by the 
shifts \n) — ► |n — 1) and \n) — ► \n + 1), respectively 

For the case of a generic initial state the algorithm used 
to generate a realization of the process may conveniently 
be expressed in term of the wave function ip(P) = (P|V>) 
in the momentum representation. In fact, for a state 
\ip(to)) = \ip) the total transition rate given by Eq. (fl~2"J) 
takes the form 



Determination of the total transition rate 



r(|^»= / dPr(P)MP) 



where 



r(p) = 



dQ<r(Q)S(Q,P) 



(19) 



(20) 



is the total rate for transitions out of a state characterized 
by the momentum P. The deterministic time evolution in 
between the quantum jumps determined by the effective 
Hamiltonian is explicitly given by 

m + r)) = J — )l , (21) 

JdPe-rO^lVKP)! 2 

while the waiting time distribution defined by Eq. (TT 
becomes 



f(t) = 1- J dP e - r ( p > T |v3(P)| 2 . 

The jumps are described by 

JdP v /^(Q,P)^(P)|P + Q> 



(22) 



(23) 



7dP5(Q,P)|^(P)| 2 
where the momentum transfers Q follow the distribution 
- fdP<r(Q)S(Q,P)\j(P)\ 2 



R(Q) 



fdPT(P)\iP(PW 



(24) 



As it immediately appears from Eqs. (|2lj) - (f24|) the choice 
of an initial state given by a finite linear superposition of 
improper momentum eigenvectors leads to very impor- 
tant simplifications (see Sec. HVBj) . 



The total transition rate (|20|) plays an important role 
in the simulation algorithm. It can be analytically calcu- 
lated in several interesting cases, e. g. considering a gas 
of free particles described by Maxwell-Boltzmann statis- 
tics so that 



r(P) = !!s» JEt f d Qa(Q) 
mi V 2tt J 

p (2mE(Q,P) + Q 2 ) 2 



(25) 



8m Q 2 
It is of great advantage to introduce the scaled momenta 



and 



K 



U 



Q 



Mv„ 



(26) 



(27) 



where v mp = yj2/mf3 is the most probable velocity of the 
gas particles. In terms of these quantities the function 
F, now expressed in terms of the scaled momentum U, 
becomes 



T(U) 



a p2V7T / 

Jo 



dK K a(K) 



d£, exp 



K 

T 



us, 



(28) 



where the variable £ denotes the cosine between the vec- 
tors U and K, and the scattering cross section a(K) has 
been supposed to depend only on the modulus of the 
momentum transfer. Integrating over £ one has 



T(U) = n eas v mp - / dKKa(K) 
U Jo 

elf 



r ! - erf I j-U 



(29) 



For the case of a constant scattering cross section, 

a(K) = a = const, (30) 
the total transition rate is found to be 



r(t/) - r I [l + 2U 2 ] 



erf(U) 
2U 



1 



where the quantity 



ro = "■ga S l'mp47rCT 



(31) 



(32) 



represents the total scattering rate corresponding to an 
incoming flux of particles with the most probable velocity 
v mp . For a Gaussian scattering cross section of the form 



a(K) = ere 



-aK 2 /4 



(33) 
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FIG. 1: (Color online) The total transition rate V(U) for a 
constant cross section [Eq. (|31[) ] (continuous line) and for a 
Gaussian cross section with o = 1 [Eq. (I34|l ] (broken line). 



one has instead 



r(U) = ?±{nUU) 
au 



The form of these functions is illustrated in Fig. Q] One 
observes that for large U the function T(U) increases lin- 
early with U in the case of a constant cross section, while 
it decreases as U^ 1 in the case of a Gaussian cross sec- 
tion. We note also that in terms of the scaled momentum 
variables the canonical mean value (P 2 /2M) cq = 3/2/3 of 
the kinetic energy of the test particle reached at thermal 
equilibrium corresponds to (U 2 ) oq — 3m/2AI . 



IV. NUMERICAL ALGORITHMS AND 
SIMULATION RESULTS 

A. Dissipative effects 

We first study relaxation to thermal equilibrium by 
considering the dynamics of an ensemble of momentum 
eigenstates of a test particle in a Maxwell-Boltzmann gas. 
Hence, using the scaled momentum variables introduced 
in Sec. [II] we investigate initial states of the form 



\m) = iu(o)>. 



(35) 



The application of the quantum jump unravelling de- 
scribed in Sec. IIIII then leads to a classical stochastic 
process U(t) for the test particle momentum. This pro- 
cess is a pure jump process the realizations of which are 
obtained through the algorithm described below. Note 
that this algorithm corresponds to the standard algo- 
rithm that is used for the stochastic simulation of classi- 
cal Markovian master equations [23[. In order to study 
relaxation to a Gaussian thermal state we will consider 
the behavior in time of first and second moments of the 
momentum distribution, as well as of various cumulants 
of the distribution. 



1. Simulation method 

The state U(0) at the initial time t = is to be drawn 
from a given distribution for the initial data. Suppose 
a jump occurred at some time t leading to the state 
TJ(to) = U. The next jump will then take place at time 
to + t, where r is a stochastic time step which is given 
by 



1 



r(co 



In 77. 



(36) 



77 is a random number which is uniformly distributed over 
the interval (0, 1), and T(U) represents the total transi- 
tion rate. Since the process is a pure jump process, U 
stays constant between to and to + T - 

At time to + r one carries out a jump by replacing 



U 



u 



AT 



K. 



(37) 



The momentum transfer K is determined as follows. 
First, one draws random numbers (K, £) that follow the 
joint probability density R{K,^) which is normalized as 



dK I d£R{K,0 = 1. 



(38) 



K is the size of the momentum transfer and £ the cosine 
of the angle between K and U. In the Monte Carlo sim- 
ulations shown below we have used the rejection method 
to determine (K, £). Second, one draws a uniformly dis- 
tributed random unit vector e. Then, the momentum 
transfer is given by the formula 



K = K 



U / - U x e . 

^ = ^- + K^e-^- (39) 



K|| is the component of K which is parallel to the particle 
momentum U, and Kj_ its component perpendicular to 
it. Repeating these steps until the desired final time t/ 
is reached one obtains a realization of the process U(t) 
over the whole time interval [0, tf\. 



2. Constant scattering cross section 

We first address momentum and energy relaxation for 
the case of a constant scattering cross section a. The cor- 
responding total transition rate T(U) is given by Eq. (|3ip . 
while the joint probability of (K, £) takes the form 



R(K,0 



To 



20?r(io 



K exp 



K 



(40) 



This probability density is illustrated in Fig.[2J For small 
U the density R(K, £) is nearly uniform in £, correspond- 
ing to an isotropic distribution, while it has a pronounced 
maximum at £ = — 1 if U is not small. In the latter 
case there is thus a strong tendency that the momentum 



FIG. 2: (Color online) The probability density R(K,£) given 
by Eq. (gOJ for U = 1. 




FIG. 3: (Color online) A single realization of the process U(t) 
for m/M = 1. 



5 10 15 20 

V 

FIG. 4: (Color online) Averages over 10 4 realizations for a 
constant cross section and m/M = 0.1. The broken lines 
represent the approximate relaxation dynamics according to 
Eqs. (|4l) and (l42|) . respectively. 




transfer K is opposite to the direction of the particle mo- 
mentum U. 

Figure [3] shows a single realization of the process U(t) 
as a 3D plot. One observes that already a few momen- 
tum kicks drive the test particle into the vicinity of the 
equilibrium value. Statistical estimates for the quantities 
(U(i)) 2 and (U(i) 2 ) for various values of m/M are shown 
in Figs. HJ 03 and [5] In these simulations we have taken a 
sharp initial state proportional to (1,0,0). We see a nice 
relaxation to the respective equilibrium values (U) 2 q = 



cq 



3m/2M for all parameter combinations 



and (U 2 ; 
used. 

A simple approximation of the relaxation dynamics can 
be obtained in the limiting case m/M <C 1. According 
to Ref. [2J| one then finds 



(V(t)) 2 « (u(o))V 



(41) 



and 



(U 2 (t)) « [<U 2 (0)) - <U 2 ) cq ] e-^ + <U 2 ) cq , (42) 
where the relaxation rate is given by 

16 m 



3^ttM 



r . 



(43) 



2 3 



FIG. 5: (Color online) The same as Fig. U for m/M = 1. 



We see from the figures that the dynamics of the mean 
squared momentum (U 2 (i)) strongly deviates from these 
approximations if m/M is not small, while the behav- 
ior of the squared mean momentum (U(i)) 2 is still well 
approximated by Eg . (|4"3"|) . A further characteristic fea- 
ture is that for m/M 3> 1 the squared mean momen- 
tum (XJ(t)) 2 decays much faster than the mean squared 
momentum (U 2 (i)). This means that the relaxation of 
the momentum and of the kinetic energy of the particle 
are characterized by different decay times, by contrast to 
typical master equations for quantum Brownian motion. 



3. Gaussian scattering cross section 

For a Gaussian scattering cross section of the form 
cr{K) = crexp[— aK 2 /4\ the total transition rate was 
given in Eq. (|34|) . The corresponding joint probability 
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FIG. 6: (Color online) The same as Fig. g] for m/M = 10. 



density becomes 



r„ 



2y/ZT(U) 



ifexp 



K 
~2 



US, 



-K 1 



(44) 

The qualitative features of this density are the same as 
for the case of a constant cross section. Simulation results 
for m/M = 1 and a — 1 are shown in Fig. [71 As expected 
we see the relaxation to the equilibrium values predicted 
by the stationary solution. In the case m/M <C 1 the 
approximations given by Eqs. (|4"Tj) and (|4"2"|) are again 
valid, where the relaxation rate now takes the form 



1R 



16 m 



M (1 + a) : 



:1V 



(45) 
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4- Cumulants of higher order 

The stochastic process U(t) describing the particle mo- 
mentum is nonlinear in the sense that large non-Gaussian 
fluctuations are dynamically generated. Starting from a 
Gaussian initial state we end up for long times with a 
Gaussian final (equilibrium) state. However, for interme- 
diate times one observes strong deviations from Gaussian 
statistics. 

To investigate these deviations one has to consider 
higher moments of the components Ui(t), i — 1,2,3, of 
the test particle momentum. To this aim we have deter- 
mined the time-dependence of the cumulants of second, 
third and fourth order (summed over the components of 
the momentum), 



K 2 



«3 



£<(Ci-<0i» 2 >, (46) 

i 

£((^-<^» 3 >, (47) 

i 

£ [((U t - (C/,)) 4 ) - 3 ((Ut (CA)) 2 ) 2 ] . (48) 



For a Gaussian distribution all cumulants of order larger 
than two vanish identically. Simulation results are shown 
in Fig. [5] We indeed see the emergence of large non- 
Gaussian fluctuations for which the cumulants K3 and K4 
are of the same order of magnitude as the variance «2- 
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FIG. 8: (Color online) Cumulants of second, third and fourth 
order obtained from averages over 10 s realizations for a con- 
stant cross section and m/M = 1. 



Decoherence effects 



FIG. 7: (Color online) Averages over 10 4 realizations for a 
Gaussian cross section with m/M = 1 and a = 1. The broken 
lines represent the approximate relaxation dynamics with the 
rate given by Eq. (|45p . 



As mentioned in Sec . HI1 diagonal matrix elements in the 
momentum representation of the quantum linear Boltz- 
mann equation coincide with the classical expression. 
Typical quantum features are therefore linked to the off- 
diagonal matrix elements and their behavior in time. The 
suppression of these matrix elements corresponds to a 
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transition from the quantum to the classical regime. It is 
therefore of interest to consider the evolution in time of a 
superposition state. As argued in Sec. IIHI superpositions 
of momentum eigenstates are preserved in the course of 
the time evolution. This implies that an initial state of 
the form 

N 

|^,(0)) =^0,(0)1^(0)), (49) 

i=l 

where the amplitudes satisfy the normalization condition 
J2iLi l a i(0)| 2 — 1j can be expressed at time t as 

N 

i^(t))=x;«i(*)iu i (t)). (so) 

z=l 

The stochastic state vector is therefore uniquely fixed 
by N momenta {Uj(t)}j = i...jv and N complex ampli- 
tudes {a.i(t)}i = \,,,N obeying the normalization condition 

£*Li l«i(*)| 2 = 1 for an y 



1. Simulation method 



where the factors fi are given by 

-±(K/2+K-Ui(t )/K) 2 

fi = , (57) 

JEj \Oj(t + T )|2 e -(K/2+K.U 3 ( t0 )/K) 2 

The momentum transfer K in these formulas is to be 
drawn from the corresponding probability density i?(K). 
To determine K one proceeds as follows. First, one draws 
an index i £ {1, 2, ... , N} with probability 

_ MMI 2 e- r ^fa)>-r(L^ )) ^ 

Given i one then draws a momentum transfer K that 
follows the probability density R(K, £,) [see Eq. (j40|) ]. 
where now £j represents the cosine of the angle between 
K and U;(t )- 

2. Decay of coherences 

We have used the simulation algorithm described 
above to study the loss of coherence of an initial state 
of the form 



Suppose that at time to a jump into the state 



|V(*0)> =5><(*o)l U i(*o)> 



(51) 



i=l 



occurred. The deterministic time evolution before the 
next jump generated by the effective Hamiltonian (jlip is 
described by 



N 



\i>(fo + °i(fo + r)|Ui(to)). (52) 



i=l 



Hence, the momenta U, stay constant while the dynamics 
of the amplitudes is given by 



a»(*o +t) = 



e -iE(Ui(t ))T /h e -T(U t (t a ))r /2 



&i(t ), (53) 



where E{U) = P 2 /2M = Mv 2 ap U 2 /2. The waiting time 
distribution F(t) for the next jump is now represented 
by a sum of exponential functions, 



N 



F(r)=l-^| ai (t ) 



2--r(£/ 4 (t ))r 



(54) 



i=l 



The jump at time t = to + r is described by the replace- 
ments 



U<(t ) — > U 4 (t ) + ^K, (55) 
ai(to + T) — > fiOti(t +T), (56) 



|^(0)) = ar^lU^O)) + a 2 (0)|U 2 (0)), (59) 

given by a superposition of two momentum eigenvectors. 
The initial state is to be drawn from a given initial distri- 
bution for the momenta Ui,2(0) and amplitudes 01,2(0). 
In the simulations shown below we have taken sharp op- 
posite initial momenta, 



U 1 (0) = -U 2 (0) = U 0) 



and equal amplitudes, 



ai(0) = 0:2(0) = 



V2' 



(60) 



(61) 



Such an initial state is a balanced coherent superposition 
of two momentum eigenvectors separated by twice Uq = 
AP/j) mp , where p mp = mv mp is the most probable mo- 
mentum of the gas particle at temperature T = 1/ksP- 
In order to study quantitatively the loss of coherence we 
have used the simulation algorithm to estimate the ex- 
pectation value 



C(t) = E 



\ ai (t)a*(t)\ 



K (OK (0)| 



2E[|ai (i)a£(i)|]. (62) 



This quantity represents the average of the absolute val- 
ues of the coefficients in front of the off-diagonal matrix 
elements of the test particle's statistical operator in the 
momentum representation, divided by its initial value. 
Hence, C(t) provides a measure for the degree of the co- 
herence of the state of the test particle. An example for 
the dynamical behavior of C(t) is shown in Fig. [5] In 
this figure we have used a sharp initial state given by 
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010 




exploiting Eq. (|40| . 



FIG. 9: (Color online) Semilogarithmic plot of the coher- 
ence C(t) obtained from an average over 10 realizations for 
m/M = 1 (constant cross section). 



Uo = (0,0,4). We clearly see an exponential decay of 
the coherence C(t) over several orders of magnitude. 

Assuming that an exponential decay of the coherence 
holds true, we define the decoherence rate jd by means 
of 



C(t) 



-lot 



(63) 



An analytical approximation for can be found with 
the help of the following argument. We take a sufficiently 
small time t such that we can approximate 

C{t) Hl-7 fl twl- T{U )t + T{U Q )t{hf 2 ). (64) 

Here, T{Ui{Q)) = T(U 2 (0)) = T(U ) represents the to- 
tal rate for a transition out of the given initial state 
([59| . Hence, T(Uo)t represents the probability for a jump 
within time t, while 1 — T(Uo)t is the probability that 
no jump occurs. Using Eqs. (f6TJ|) and (|6"Tj) we see from 
Eq. ([53]) that C does not change during the deterministic 
drift. On the other hand, if a jump does occur then C 
changes from its initial value C = 1 to C = /1/2 as may 
be seen from Eq. (p)6"|) . According to Eq. (|57p we have 



hh = 



(65) 



Thus, Eq. ([64]) represents the change of C as a result 
of two alternatives, namely that a jump does occur or 
that it does not (for small enough t we can have at most 
one jump). Finally, (/1/2) denotes the average of /1/2 
taken over the possible momentum transfers during the 
first jump. We therefore get 



7 d « r([/ )(l - (A/a)) 

= r(C/ )(l-Sech(K-U )) 



(66) 



This formula can be analytically evaluated in several 
cases. For a constant scattering cross section one finds 



Id = T(U ) - T 



erf(Efr) 
Uo ' 



(67) 



Fig. [TU] demonstrates the extremely good agreement be- 
tween the stochastic simulation results and this analytical 
approximation. The Monte Carlo estimates for the deco- 
herence rate have been obtained by least squares fits to 
the simulation data within the time interval in which the 
coherence decays to 1% of its initial value. For all cases 
shown we find a nearly perfect exponential decay. Our 
argument shows that it is just the real scattering events 
that mainly cause the decoherence during the early phase 
of the dynamics, by contrast to the virtual transitions de- 
scribed by the non-Hcrmitian drift Hamiltonian. 




FIG. 10: (Color online) The decoherence rate 70 in units of 
To as a function of the initial momentum Uo- Points: Least 
squares fits of the simulation data for m/M = 1. Continuous 
line: Analytical estimate given by Eq. (1671 1, 



0.25 




FIG. 11: (Color online) The same as Fig. [10] for a Gaussian 
cross section with a = 1. The continuous line represents the 
analytical estimate given by Eq. 



For a Gaussian cross section Eq. ([6"6")) leads to 

erf(£/ ) 1 



ID = r(u a ) - To 



Uo 1 + a ' 



(68) 
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In Fig. [TT] we compare this expression with the simula- 
tion results. Although the quantitative agreement is ob- 
viously not as good as in the case of a constant cross sec- 
tion, the analytical formula (j6"8)) still yields a reasonable 
estimate for the order of magnitude of the decoherence 
rate. 

The above results enable us to study the relationship 
between the time scales characterizing the different phys- 
ical phenomena. For the case oim/M <C 1 the time scale 
for relaxation is given by Eq. (|43[) . Considering a wide 
separation in momentum of the initial superposition on 
the scale set by the momentum of the gas particles, so 
that Uo — AP/p mp ^> 1, one can consider the corre- 
sponding limiting value TqIIo of Eq. (|67|) . and the ratio 
of the decoherence rate to the relaxation rate is given by 



Yd 
Yr 



VttM 
16 m 



Uo » 1. 



(69) 



This shows that the decoherence time j^ 1 for the loss of 
coherence in momentum space can be much smaller than 
the relaxation time Yrj i- e. than the time it takes for 
the relaxation of the energy of the test particle. 

It is also of interest to compare the relaxation timescale 
with the decoherence rate tjd of a superposition of po- 
sition eigenstates, widely separated on the scale set by 
the thermal wavelength of the gas particles, which corre- 
sponds to H divided by the typical momentum transfer, 
so that AX/Ath ~> 1- In such a case for m/M -C 1 
the decoherence rate is set by the total transition rate 
[E0, HE HE USUI] , so that for a test particle much slower 
than the gas particles one has r\u w 2To/y/ir. Thus, the 
ratio 



UR 

IR 



3M 



> 1 



(70) 



shows again that the decoherence time for the loss of 
coherence in position space can be much smaller than the 
relaxation time Yn > but still longer than the decoherence 
time in momentum space j^ 1 . 



V. CONCLUSIONS 

We have developed a stochastic unravelling of the 
quantum linear Boltzmann equation which leads to an 
efficient Monte Carlo simulation technique, despite the 
appearance in the equation itself of quite complicated 
operator-valued expressions. The latter are responsible 
for deviations from Gaussian statistics, at variance with 
typical master equations used for the description of quan- 
tum Brownian motion. A crucial feature of the method 
is that the developed algorithms fully exploit the trans- 
lation covariance and thus allow full three-dimensional 
stochastic simulations of the quantum Boltzmann equa- 
tion. In particular, the method does not require the 



introduction of a discretization in momentum space, so 
that the continuous sum over the Lindblad operators in 
Eq. ([9]) can be exactly accounted for. 

The method proposed here suggests many further 
physically relevant applications. For example, one can 
extend the algorithm to the regime where effects from 
the Bose or Fermi statistics of the quantum gas come 
into play. In fact, starting from Eqs. ((3]) and {5| one 
can analytically work out the corresponding expressions 
of the dynamic structure factor for a free gas of particles 
obeying Bose or Fermi statistics, coming to [4[ 



S b /f(Q,E) — 



1 



2nm 2 



(2^)3 n gas f3Q 1 - e? E 



x In 



1 =p z exp 


r 


(2mE+Q 2 ) 2 ' 




8m 


Q 2 




1 =p z exp 





(2mE-Q 2 ) 2 ' 




8m 


Q 2 





where the upper signs refer to the Bose case and the lower 
signs to the Fermi case, and z denotes the fugacity of the 
gas. By use of this expression the algorithm thus enables 
Monte Carlo simulations of the behavior of test particles 
in a Bose or Fermi gas to identify genuine effects of the 
quantum statistics. More generally, the method may also 
be applied to an interacting quantum gas provided an (at 
least approximate) expression for the dynamic structure 
factor 5(Q, P) is known. 

A further example is the investigation of the important 
problem of decoherence in position space. This can be 
done by use of initial states representing superpositions 
of localized wave packets, with the aim of determining 
the corresponding position space decoherence time scales. 
Localized wave packets may of course be represented by 
introducing a discretization of position space. However, 
it seems that it is much more efficient to invoke the trans- 
lation covariance and to describe spatially localized states 
by appropriate superpositions of momentum eigenstates 
as discussed in Sec. IIVB1 or, more generally, by super- 
positions of wave packets localized in momentum space. 
Further examples of application include the extension of 
the method to the determination of multitime correlation 
functions, to the treatment of particles with internal de- 
grees of freedom, and to the case of an operator-valued 
scattering amplitude, i. e. to the case that the scattering 
cross section depends on the momentum of the incoming 
test particle. 
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